MULTILEVEL SPARSE KERNEL-BASED INTERPOLATION 
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Abstract. A multilevel kernel-based interpolation method, suitable for moderately high-dim- 
ensional function interpolation problems, is proposed. The method, termed multilevel sparse kernel- 
based interpolation (MLSKI, for short), uses both level-wise and direction-wise multilevel decomposi- 
tion of structured (or mildly unstructured) interpolation data sites in conjunction with the application 
of kernel-based interpolants with different scaling in each direction. The multilevel interpolation al- 
gorithm is based on a hierarchical decomposition of the data sites, whereby at each level the detail is 
added to the interpolant by interpolating the resulting residual of the previous level. On each level, 
anisotropic radial basis functions are used for solving a number of small interpolation problems, 
which are subsequently linearly combined to produce the interpolant. MLSKI can be viewed as an 
extension of d-boolean interpolation (which is closely related to ideas in sparse grid and hyperbolic 
crosses literature) to kernel-based functions, within the hierarchical multilevel framework to achieve 
accelerated convergence. Numerical experiments suggest that the new algorithm is numerically stable 
and efficient for the reconstruction of large data in R d X M, for d = 2, 3, 4, with tens or even hundreds 
of thousands data points. Also, MLSKI appears to be generally superior over classical radial basis 
function methods in terms of complexity, run time and convergence at least for large data sets. 
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1. Introduction. Over the last four decades, radial basis functions (RBFs) have 
been successfully applied to (scattered) data interpolation/approximation in t d x I 
(see, e.g., [35] and the references therein for a literature review). The interest on 
kernel-based and, in particular, on RBF interpolants can be traced in their abil- 
ity to produce global interpolants of user-defined smoothness without the shortcom- 
ings of multivariate polynomial interpolation. These interpolants admit generally 
good convergence properties and they can be implemented in (essentially) dimension- 
independent fashion, making them potentially attractive for a number of applications. 

Despite the above attractive properties, RBF interpolation can be cumbersome 
in practice. Solving the resulting linear system is challenging due to both the density 
and the ill-conditioning of the resulting interpolation matrix (see, e.g., [121 E2])- A 
number of techniques have been proposed to deal with the ill-conditioning of the 
interpolation system. 

For smooth basis functions such as the Gaussian and the multiquadric, the ill- 
conditioning depends crucially on a width parameter. Novel QR-based algorithms 
which eliminate the ill-conditioning problem for thousands of points have been devel- 
oped [201 [T51 H3] . An alternative family of methods, where small local problems are 
solved which generate approximate cardinal functions has also been proposed [THl [3] . 
The matrix associated with such an approach is much better conditioned, and al- 
lows for rapid solution using iterative methods. Finally, in (40] . a stable interpolation 
algorithm on carefully selected nodes for Gaussian basis functions is given. 

The introduction of RBFs with compact support [501 HZ] aims to address the 
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density issue of the interpolation matrix. Moreover, a number of techniques have 
been developed to reduce the complexity of calculating the interpolant, involving 
multipole type expansion for a variety of RBFs [3J . Using such methods is possible to 
compute an RBF interpolant with 0(k log fc)-computations for quasi-uniform data, 
where k is the number of data sites though, to the best of our knowledge, we are not 
aware of any methods which guarantee a bounded number of iterations independently 
of k. 

Thus, the complexity of kernel-based interpolation remains a challenge when d > 
2, as the number of data points is required to grow exponentially with respect to the 
dimension d to ensure good convergence rates. It is not surprising, therefore, that the 
use of RBFs in practice has been largely limited to tens of thousands of data sites, 
which, in turn, has restricted their application to low d, typically d = 1,2 or 3. 

This work is concerned with the introduction of a kernel-based interpolation 
method on structured or mildly unstructured data sites, which aims to address the 
computational complexity issues of RBF interpolation for d > 2, while simultaneously 
reducing the ill-conditioning of the resulting interpolation problem. The new scheme, 
termed multilevel sparse kernel-based interpolation (MLSKI), is based on a hierar- 
chical decomposition of the data sites, whereby at each level the detail is added to 
the interpolant by interpolating the resulting residual of the previous level. On each 
level, anisotropic radial basis functions are used for solving a number of small interpo- 
lation problems, which are subsequently linearly combined to produce the interpolant; 
the new method can be viewed as an extension of (and, indeed, it has been inspired 
from) the idea of d-boolean interpolation [TOl HH ESI I3H1 I23J to kernel-based functions, 
which, in turn, is closely related to ideas in sparse grid [5TJ [S7J SI 1311 HH| and hy- 
perbolic crosses [351 [J literature. We note that in j33J hyperbolic cross products of 
one dimensional RBFs have been considered. The hierarchical multilevel framework 
used to achieve accelerated convergence is relatively standard in the RBF literature 

[IB EES EH ESI- 

In the simplest setting, the MLSKI algorithm assumes that the data sites are 
on a Cartesian uniform grid of size N d in E d , with N = 2™ + 1, for a final level 
n € N. For each level < I < n of the MLSKI algorithm, we construct a sparse 
kernel-based interpolant of the interpolation residual as follows. We consider 0(Z d_1 ) 
carefully chosen subsets of the data points of level I, each subset having size 0{2 l ) 
data points. On each of these subsets, which we shall refer to as partial grids, we 
solve the interpolation problem. As the partial grids are anisotropic in nature, we 
employ appropriate anisotropically scaled kernels (anisotropic RBFs 8, 9, 2j) for each 
interpolation problem on each partial grid. Once all the 0(^ d_1 ) interpolants on the 
partial grids have been computed they are linearly combined to give the total inter- 
polant (of the residual) for level I. Hence, the complexity of the resulting MLSKI 
algorithm is dominated by the complexity of the last step, i.e., one needs to solve 
0(n <J ^ 1 ) interpolation problems of size N and linearly combine the resulting inter- 
polants on the partial grids. This, in turn, implies that the MLSKI algorithm admits 
(at least theoretically) O (a [N) log'' -1 a(N) )-complexity, where TV is the number of 
grid points in each direction, and cr(iV) is the complexity of the univariate TV-point 
algorithm. 

The computation of each interpolant on the partial grids is completely indepen- 
dent from the other interpolants on each level I. Therefore, the MLSKI algorithm 
is ideally suited for implementation on parallel computers. The MLSKI algorithm is 
tested in practice for a number of relevant test cases with d = 2,3,4. The numerical 
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Fig. 2.1: Gaussian and anisotropic RBFs, with shape parameter c = 0.5. 



experiments suggest that MLSKI is numerically stable in practice and efficient for the 
reconstruction of large data in R d x K, for d = 2, 3, 4, with hundreds of thousands of 
data points. Also, MLSKI appears to be generally superior over classical radial basis 
function methods in terms of complexity, run time and convergence, at least for large 
data sets. 

The remaining of this paper is organized as follows. In Section [2] we discuss 
anisotropic versions of radial basis functions, suitable for interpolation on data sites 
with anisotropic characteristics. In Section [3j we introduce the sparse kernel-based 
interpolation method, which will be used in Section [4] at each level of the multilevel 
algorithm. The stability and implementation of the MLSKI algorithm is discussed in 
Section 
Section 



51 A series of numerical experiments is given in Section [6j for d = 2, 3, 4. In 
7 we draw some conclusions on the new developments presented in this work 
and discuss possible extensions. 

2. Anisotropic RBF interpolation. Radial basis functions are radially sym- 
metric by construction having hyper-spheres as their level surfaces. However, inter- 
polation of data with anisotropic distribution of data sites in the domain requires 
special consideration. To this end, anisotropic radial basis functions (ARBFs) have 
been introduced and used in practice [H [9j [2] ; they are also known as elliptic basis 
functions as they have hyper-ellipsoidal level surfaces. 

Definition 2.1. Let <p(\\ • — Xj||) be a given RBF centred at x^- G R d and let 
A £ R dxd be an invertible matrix. The anisotropic radial basis function ipA is defined 
bv<p A Q\.-x j \\) = <p(\\A(--jL S )\\). 

Evidently, ipA (IMI) = <p{ || • || ) when A is the d x d-identity matrix. In Figure pO] the 
Gaussian RBF and the corresponding anisotropic Gaussian RBF centred at (0.5, 0.5) 
for A = diag(2 5 ,2 2 ) are drawn. (Here, we use the notation diag(w) € M. dxd for 
a diagonal matrix with diagonal entries are given by the components of the vector 
v € R d .) 

We now discuss the the solution of an anisotropic interpolation problem. We 
restrict the discussion to positive definite kernels, which will be used later, though 
anisotropic versions of conditionally positive definite REBs is completely completely 
analogous. 

For data sites X := {xi, . . . ,xjv} contained in a computational domain ft C K d , 
we consider the interpolation data {(xj,j/j) : £ Xi = 1, ...,iV}. Let ip be a 
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positive definite radial function and A £ R dxd be a chosen invertible matrix. Then, 
the anisotropic RBF interpolant 5a is defined by: 

N 

^(x) = ^c^ A (||x-x J -||) ) xGtt, (2.1) 

i=i 

for Ci such that the interpolation conditions 5a(x^) = t/j, i = 1, . . . ,N, are satisfied. 

The well-posedness of the interpolation problem is guaranteed (for positive defi- 
nite kernels) as a direct consequence of the invertibility of the scaling matrix A. We 
refer to [2] for the error analysis of anisotropic RBF interpolation. 

3. Sparse kernel-based interpolation. We describe the basic sparse kernel- 
based interpolation (SKI) method that will be used in each step of the multilevel 
algorithm. The SKI method can be used also as a single step interpolation method. 

The starting point is the observation that, assuming sufficient smoothness of the 
interpolation data, the number of points required to provide a given accuracy can be 
dramatically reduced when basis functions with carefully constructed direction-wise 
anisotropic scaling are used. This way, notwithstanding the approximation strength 
coming from a few directions only, there is only negligible loss of accuracy, due to the 
additional smoothness assumed. Hyperbolic cross products and sparse grids use this 
idea in the context of high dimensional numerical quadrature, approximation and in 
numerical solution of partial differential equations. 

In [13] , hyperbolic crosses of tensor products of one-dimensional RBFs have been 
considered without numerical assessment of the resulting method. The use of non- 
tensor product d-dimensional RBFs in the hyperbolic-cross/sparse-grid setting is not 
straightforward, as such approximation spaces are characterised by basis functions 
with differently anistropic scalings in different directions. Using such scalings would 
not guarantee the well-posedness of the resulting kernel-based interpolation problem. 

An alternative approach can be motivated by the, so-called, d-boolean interpo- 
lation [TU] or, the related combination technique [27J H2, for piecewise polynomial 
interpolation on sparse grids. A key point in this approach is that the piecewise poly- 
nomial interpolant is equivalent to a linear combination of interpolants constructed 
on a carefully selected collection of subsets of the (anisotropic) basis functions. Each 
such subset consists of translations of the same (anisotropic) basis function. 

To construct the sparse kernel-based interpolant, we solve a number of anisotropic 
radial basis function interpolation problems on appropriately selected sub-grids and 
we linearly combine the resulting partial interpolants to obtain the sparse kernel-based 
interpolant. 

More specifically, let := [0, l] d , and let u : £1 — > R. The extension to general 
axiparallel domains is straightforward; we refrain from doing so here in the interest of 
simplicity of the presentation only. Comments on possible extensions to more general 
domains will be given in Section [7J 

For a multi-index 1 = (^,...,^) g N d , we define the family of directionally 
uniform grids {Xi : 1 G N d }, in Q, with meshsize h\ = 2" 1 := (2~ h , . . . ,2~ ld ). That 
is, Xi consists of the points xu := (x^^ , . . . , xi d ,i d ), with xi^. = ij2 3 , for ij — 
0,1, . . . ,2 lj , j = 1, . . . ,d. The number of nodes TV 1 in Xi is given by 

d 

N 1 = l[{2 1 + 1). 

i=l 
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Fig. 3.1: Sparse grid X 4 ' 2 via (3.1 1. 



If hi i —2 ", for all i = 1, • • • , d, Xi is the uniform full grid of level n, having size 
N=(2 n + l) d ; this will be denoted by X n - d . 

We also consider the following subset of X n ' d , 

X n ' d := (J Xi, (3.1) 

|i| 1=n+ (d_i) 

with := li + • • • + la, which will be referred to as the sparse grid of level n in d 
dimensions. We refer to Figure 3.1 for a visual representation of (3.1) for n = 4 and 
d — 2. Notice that there is some redundancy in this definition as some grid points are 
included in more than one sub-grid. 

We want to evaluate the interpolant at the constituent sub-grids Xi. As the 
constituent grids admit different density in each coordinate direction, we shall make 
use the anisotropic RBFs from Section [2] To this end, for each multi-index 1 = 
(li, . . . , Id), we define the transformation matrix A\ £ M. dxd by 

A, :=diag(2'V..,2^). 

The anisotropic RBF interpolant S^, of u at the points of Xi is then defined by 

N 1 

S Al (x):=Y, c M\\Mx-^)\\)> (3-2) 
i=i 

for x £ f], where Cj £ E are chosen so that the interpolation conditions 
are satisfied. 

To construct the sparse kernel-based interpolant (SKI, for short) on the sparse 
grid X™' d , the sub-grid interpolants Sa { are linearly combined using the formula 

s n (x) = 5j(-i)»( d ~ 1 ) E ( 3 - 3 ) 



The combination formula (3.3) has been used in the context of d-boolean Lagrange 
polynomial interpolation |10j. and in the combination technique for the numerical 
solution of elliptic partial differential equations using the finite element method on 
sparse grids [371 121] ■ For d — 2, (|3.3[) becomes 



|l|i=n+l |l|i=n 



(3.4) 
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Fig. 3.2: The construction of the sparse kernel-based S4 interpolant on X 4 ' 2 . 



that is, the first term on the right-hand side of (3.4) gives the sum of interpolants 
on the sub-grids of level n + 1, while the second term on the right-hand side of (3.4) 
subtracts the redundant points visited more than once. We refer to Figure [3~2] for an 
illustration when d — 2 and n — 4. 

Hence, SKI uses a dimension-wise multilevel decomposition of interpolation data 
sites in conjunction with the application of kernel-based interpolants Sai(-) with dif- 
ferent scaling in each direction. The sparse kernel-based interpolant S n can be imple- 
mented in a quite straightforward fashion by utilising existing, fast, RBF interpolation 
algorithms: the only modification needed is the introduction of a scaling for each sub- 
grid problem. We note that each interpolation problem can be solved completely 
independently, rendering the resulting SKI method ideally suited for implementation 
in parallel computers. 

Moreover, we observe the size of each sub-grid problem is 0(2"), where n is the 
number of levels, i.e., it is independent of the dimension d. One has to solve 0{n d ^ 1 ) 
such sub-grid problems to obtain the sparse kernel-based interpolant. Observe that 
0(2") is the number of points of the corresponding full grid X" in each space direction. 
Hence, if the rate of convergence for the SKI algorithm is comparable with the one 
of the standard RBF algorithm for a sufficiently large class of underlying functions 
u, the benefits in complexity are potentially significant, especially for d > 2. Indeed, 
when the underlying function u is assumed to admit sufficient regularity (e.g., if the 
regularity of u is characterised by anisotropic Sobolev spaces) [lOl [27l [22] and [43] the 
rate of convergence on sparse grids for tensor-product of piece-wise polynomials and 
of one-dimensional RBFs, respectively, is shown to be optimal (modulo a logarithmic 
factor) . Although, there is no general proof at this point that this is also the case for 
the SKI method presented here, numerical experiments presented below show that 
for the multilevel version of the SKI algorithm, described in the next section, good 
convergence results are also observed. We finally remark that for the case of Gaussian 
interpolation, due its tensor-product nature, the theoretical results of |43j on tensor- 
product one-dimensional sparse RBF interpolation should be applicable to the case 
of SKI also. 

4. Multilevel Sparse Kernel-Based Interpolation. Multilevel methods for 
RBFs [171 [331 113 HH [HI EH, ESI [SI El] combine the advantages of stationary and 
non-stationary interpolation, aiming to accelerate convergence and to improve the nu- 
merical stability of the interpolation procedure. The basic idea of multilevel methods 
in this context is to interpolate the data at the coarsest level, and then update by 
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Fig. 4.1: Nested-ness of sparse grids X™' 2 , n = 1, 2, . . . , 6. 



interpolating the residuals on progressively finer data sets using appropriately scaled 
basis functions. 

The setting of SKI is naturally suited to be used within a multilevel interpolation 
algorithm. Indeed, the sparse grids fr om l ower to higher level are nested, i.e., X n,d C 

for an illustration when d = 2. Moreover, 



4.1 



r n +i,d £ or n g pj. we re f er to Figure 
each sub-grid interpolant uses appropriately scaled anisotropic basis function with 
the scaling being proportional to density of the corresponding constituent sub-grid. 
Finally, due to the geometrical progression in the problem size from one sparse grid to 
the next, a multilevel algorithm would not affect adversely the attractive complexity 
properties of SKI. 

The multilevel SKI (MLSKI, for short) algorithm is initialised by computing the 
SKI S n(j at the coarsest designated sparse grid X n °- d and set A := S n . Then, 
for k = 1, ...n, we compute to be the sparse grid interpolant of the residual 
u — YljZo A? 011 ^- k ' d - The resulting multilevel sparse kernel based interpolant is then 
given by 

n 
3=0 



5. Stability of SKI and MLSKI. The numerical stability of RBF interpola- 
tion is a challenging issue, due to the ill-conditioning of the respective interpolation 
matrices in standard bases; we refer to [2H1 H] and the references therein for a dis- 
cussion. It is, therefore, evident that the numerical stability of SKI and MLSKI 
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algorithms can be assessed by considering the maximum condition number of the 
constituent sub-grid interpolation problems. 

The main factors affecting the conditioning of an RBF interpolation problem are 
the problem size, the separation distance qx of the respective data set X, given by 

q x := ^min||xi - Xj-|| 2 , Xj.Xjgl, 

Z 1,3 

and the shape parameters of the RBF (which adjusts the "strength" or the support, 
depending on the type of RBF). The choice of shape parameter and the monitoring of 
the separation distance require more attention in practice than the problem size, which 
is of secondary importance [49j[T2]. Recipes for choosing RBF shape parameters in 
various contexts have been presented, e.g., in [25J [5TJ [5J [HI [551 [7] , while the dependence 
of the stability on the shape has been addressed in (32 [201 HH1 HH1 H3J HE] , where stable 
algorithms for RBF interpolation with small shape parameters have been proposed. 

SKI and MLSKI naturally include scaling adjustment at each level. So, the con- 
dition number of each sub-grid interpolation matrix appears to be (mildly) affected 
by the sub-grid size, as we shall see in Section [6] This is, indeed, a very attractive 
feature of the SKI and MLSKI algorithms as the growth of degrees of freedom on each 
sub-grid grows independently of the problem dimension d. 

Observing that, due to the anisotropy of the scaling matrices A\, the separation 
distance on each pulled-back (multiplied by A^ 1 ) sub-grid is constant with respect to 
the level n. This suggests that we can use standard heuristics for the choice of shape- 
parameters for multilevel RBF interpolation (see (17 ] 155 ] 158 1 11 ^ PI 158 1 1 ^ 154 ^ 155] . 
or |49) for a review). In particular, we set 

c =l? ' t 5 - 1 ) 

Kqv- n ,d 



for some constant K > 0. Extensive numerical experiments presented in [46j suggest 
that the choice K = 1, gives very small condition numbers at the cost of (mildly) lower 
accuracy of the interpolation problem. The choice K — 3, produces larger but mostly 
safe condition numbers (in this work we shall refer to a condition number k as being 
safe when n < 10 10 ), with the resulting computations admitting good convergence 
properties. For K = 3, we have 0.2 < c < 0.8. For K > 3, the ill-conditioning 
gradually increases with K, resulting in unstable computations. 

6. Numerical Experiments. We present a collection of numerical experiments 
for d = 2,3,4, where the implementation of SKI and MLSKI algorithms is assessed 
and compared against both the standard RBF interpolation on uniform full grids and 
its standard multilevel version (MLRBF); see, e.g, [35J for a review. 

The algorithm implementation has been performed in MATLAB© on a quad- core 
2.67GHz Intel Xeon X5550 CPU with 12GB of RAM, without taking advantage of the 
possibility of parallel implementation of SKI or MLSKI. The SKI/MLSKI algorithm is 
able to solve d-variate interpolation problems on sparse grids up to 114, 690 centers for 
d = 3, and up to 331, 780 centers for d — 4, respectively. The same basic interpolation 
solver for classical RBF on the same machine could only solve problems of size nearly 
15, 000 regardless of the dimension d, due to the size of the resulting linear system. 

We stress that no attempt has been made to use fast algorithms, e.g., the fast 
Gauss transform of Greengard and Strain [53] for the RBF interpolation problems 
for either the standard RBF, MLRBF or the SKI, MLSKI algorithms. The use of 
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SGnode 


Max-error 


RMS-error 


Cond. no 


Time 


9 


6.2215e-l 


1.8363e-l 


2.6912c+3 


<1 


21 


3.3237e-l 


7.6547e-2 


2.5325e+4 


<1 


49 


1.1130e-l 


3.8660e-2 


2.8184e+5 


<1 


113 


4.0379e-2 


1.0835e-2 


2.6522e+6 


<1 


257 


1.2649e-2 


2.5117e-3 


2.9516e+7 


<1 


577 


2.4678e-3 


4.0273e-4 


1.7591c+8 


<1 


1281 


2.2043e-4 


2.1030e-5 


1.0484e+9 


1 


2817 


3.5287e-5 


2.5391e-6 


2.3229e+9 


5 


6145 


6.2139e-6 


3.2696e-7 


5.1468c+9 


33 


13313 


1.1784e-6 


4.2920e-8 


6.5016c+9 


281 


28673 


2.1204e-7 


5.6557e-9 


8.2129c+9 


2397 


61441 


4.1321e-8 


7.6854e-10 


8.7056c+9 


21586 



Table 6.1: MLSKI for up2D, d = 2, using Gaussians with c = 0.45. Error evaluated 
at 25, 600 Halton points. 



fast or stable algorithms for RBF interpolation could be also used within the sub- 
grid interpolation problems of SKI/MLSKI. Indeed, any modern fast evaluation or 
numerical stabilization technique in RBF interpolation of a single problem can be 
incorporated into the SKI/MLSKI framework, possibly substantially improving the 
results presented here. 

6.1. Experiment 1. The MLSKI algorithm with Gaussian basis functions is 
applied to the (standard) benchmark Frankc's function up 2 D for d — 2, given by 

U F2D (X 1 ,X 2 ) := ^ e (-(9-i-2) 2 -(to2-2) 2 + 3 e _ ((tol + 1)2)/ 49-((9x 2 + l)^)/10 

4 , 4 , (6-1) 

+ 2 5 



In Table 6.1 "SGnode" stands for the number of sparse grid centers used, "Max- 
error" and "RMS-error' are the maximum norm and root-mean-square (L 2 -norm) 
errors, respectively, evaluated at 25,600 Halton points (37], "Cond. no" stands for 
the largest 2-norm condition number of the sub-grid interpolation matrices. Finally, 
"Time" is the time in seconds required to solve the interpolation problem in the above 
computer. We note that the tic ; and toe ; commands of MATLAB© have been used 
to produce the run times, which is accurate for longer run times. 

In Table [6~2| the MLSKI algorithm with Gaussian basis functions is applied to a 



three-dimensional version of Franke's function Up3D, given by 



™(^l,X 2 ,X 3 )=^(-^- 2 ) 2 -^- 2 ) 2 -^- 2 ) 2 )/ 4 



_1_ 3 e _((9 a;i + l) 2 )/49-((9a :2 + l) 2 )/10-((9x3 + l) 2 )/29 

\ (6-2) 

+ 2 

_ I e -((9*l-4) 2 )/4-(9:r 2 -7) 2 -((9:r3-5) 2 )_ 
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SGnodc 


Max- error 


RMS-error 


Cond. no 


Time 


27 


6.8808c-l 


1.0179e-l 


1.4863c+4 


<1 


81 


5.5853c-l 


7.7339e-2 


2.5376c+5 


<1 


225 


2.4324e-l 


3.8389e-2 


6.2956c+4 


<1 


593 


1.5884e-l 


2.1676e-2 


1.0617e+7 


8 


1505 


6.2918e-2 


6.7591e-3 


6.9125e+5 


9 


3713 


1.34036-2 


1.7755e-3 


4.2986c+8 


12 


8961 


2.2041e-3 


2.2448e-4 


6.2618c+6 


42 


21249 


3.3081e-4 


2.9755e-5 


9.2544e+9 


328 


49665 


8.9456e-5 


4.5151e-6 


2.1038c+7 


2900 


114689 


1.5829e-5 


5.7471e-7 


1.0594c+ll 


25160 



Table 6.2: MLSKI for UF3D, using Gaussians with c as in (5.1) for K 
evaluated at 125, 000 Halton points. 



3. Error 



SGnodc 


Max-error 


RMS-error 


Cond. no 


Time 


81 


7.9105e-2 


4.4589e-2 


3.6544e+5 


<1 


297 


2.4067e-2 


1.0677e-2 


8.6224e+6 


<1 


945 


1.9844c-2 


6.3598e-3 


1.0568c+6 


<1 


2769 


5.6653e-3 


1.2672e-3 


3.6076e+8 


2 


7681 


4.7096e-3 


8.2613e-4 


1.4065e+7 


14 


20481 


1.3155e-3 


1.5425e-4 


1.4848e+10 


164 


52993 


1.1548e-3 


1.0690e-4 


1.2741e+8 


1766 


133889 


3.2099e-4 


1.9243e-5 


6.0115e+ll 


16442 


331777 


2.8385e-4 


1.3934e-5 


1.1542c+9 


169643 



Table 6.3: SKI for u q . 
at 194, 481 Halton points. 



using Gaussians with c as in (5.1 1 for K = 3. Error evaluated 



Next, we consider the interpolation problem of the simple function u qua( i : [0,1] 
with 



id{xx,X2, Xz,Xi) := 4 4 J^iz; l (l - x{). 



i=l 

We employ the SKI algorithm for d = 4; the results are given in Table [673] 

6.2. Experiment 2. We now turn our attention to the conditioning of the sub- 
grid interpolation problems. In the previous experiment a, mostly safe, maximum 
condition number for each level's sub-grid interpolation problems is observed, using 
K = 3 in (5.1). Here, we investigate the choice K = 1 which gives O(l) condition 



numbers at the cost of a slight reduction in the quality of the approximation in the 
MLSKI algorithm. 

The results from the MLSKI algorithm using Gaussians with shape parameter c 



chosen from (5.1 1 with K — 1 are given in Table 6.4 The "Time" -column is ommittcd 



for brevity as it is almost identical to the one of Table |6.2| We observe that the error 
is roughly 5 times worse, but the condition numbers of the resulting sub-grid problems 
are nearly optimal. 
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SGnode 


Max-error 


RMS-error 


Cond. no 


27 


7.0968e-l 


1.0531e-l 


1.8 


81 


5.5864e-l 


7.6444e-2 


2.6 


225 


3.2513e-l 


4.8817e-2 


1.5 


593 


1.3272e-l 


1.5118e-2 


3.0 


1505 


7.8689e-2 


7.9827e-3 


1.5 


3713 


2.1970e-2 


2.1392e-3 


3.2 


8961 


1.0543e-2 


9.9965e-4 


1.6 


21249 


1.7569e-3 


1.7839e-4 


3.4 



Table 6.4: MLSKI for UF3D, using Gaussians with c as in (5.1) for K = 1. 
evaluated at 125, 000 Halton points. 



Error 



6.3. Experiment 3. We continue by comparing the SKI and MLSKI algorithms 
with standard RBF interpolation on full grids and with its (standard) multilevel 
version on full grids, henceforth denoted by MLRBF. 

In Figure [6T] the root mean-square error of RBF, MLRBF, SKI and MLSKI, 
respectively, for up3Di are plotted against the number of data sites N and against 
the computation time (in seconds). For RBF and MLRBF, N denotes the number of 
data sites on the full grid X™, whereas for SKI and MLSKI, N refers to the number 
of sparse grid nodes (i.e., N = SGnode). (We note in passing that using standard 
"isotropic" RBF interpolation on sparse grids does not appear to result to a convergent 



algorithm.) The choice of the shape parameter is given by (5.1) with K = 3 for 
SKI/MLSKI for all experiments; the shape parameter for RBF/MLRBF is chosen so 
as to have comparable condition numbers in all cases. We observe that RBF/MLRBF 
interpolation appears to perform better than SKI/MLSKI when the error is plotted 
against N. The SKI/MLSKI algorithms are, on the other hand, able to calculate 
larger problems and they do so efficiently in terms of complexity, at least for the case 
of larger problems. This is manifested in the RMS-error versus computation time plot 
in Figure |6.1| 

As discussed in Section [3] ti-Boolean Lagrange interpolation using polynomials 
requires additional smoothness of mixed derivatives in order to ensure the essentially 
optimal convergence rate. To test the extend to which this is also the case for the 
MLSKI algorithm, we consider the function ur^d : [0, l] 3 — > R, with 



UR3D 



{xi,x 2 , x 3 ) := (r 2 + r 4 ) logr, with r = {x\ + x\ + x\) 1/2 . 



This function does not posses smooth mixed derivatives of (arbitrarily) high order. 
In Figure |6.2| the root mean-square error of RBF, MLRBF, SKI and MLSKI, re- 
spectively, for urzd, are plotted against the number of data sites N and against the 
computation time (in seconds). Somewhat surprisingly, we observe that the ML- 
SKI algorithm outperforms RBF/MLRBF interpolation both in terms of degrees of 
freedom and in terms of computational time, at least for larger problems. 

Next, we turn our attention to the problem of interpolating five-dimensional data, 
i.e., d = 4. This is a challenge in practice with known methods: standard RBF 
interpolation is limited to few tens of thousands of data points. This, in turn, results 
to using O(10) data points in each coordinate direction, thereby, limiting the resulting 
approximation quality. 

We compare the four interpolation methods on a four-variate version of the 
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Fig. 6.1: Convergence of Gaussians for UF3D- Error evaluated at 125,000 Halton 
points. 
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Fig. 6.2: Convergence of Gaussians for Ubzd- Error evaluated at 125,000 Halton 
points. 



Franke's function ufad '■ [0, l] 4 — > M, with 



U F 4d(x X , ...,X A ):= \ e (-(9- 1 -2) 2 -(9^-2) 2 -(9x3-2)^)/4-(9 a;4 -2)^)/8 



4 
1 

-e 
2 

1 

-e 



-((9x 1 +l) 2 )/49-((9a; 2 + l) 2 )/10-((9a:3 + l) 2 )/29-((9a:4 + l) 2 )/39 



-((9zi-7) 2 )/4-(9z 2 -3)M(9z3-5) 2 )/2-((9z 4 -5) 2 )/4 



-((9xi-4) 2 )/4-(9s 2 -7) 2 -((9z 3 -5) 2 )-((9x 4 -5) 2 ) 
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N Time 
Fig. 6.3: Convergence of Gaussians for ufad- Error evaluated at 194,481 Halton 
points. 




10 2 10 3 10* 10 B 10 3 10 1 10 1 10 3 



N Time 
Fig. 6.4: Convergence of Gaussians for Ubad- Error evaluated at 194,481 Halton 
points. 



as well as the four-variate version ur/±d of ur3d, which is defined completely analo- 
gously. The convergence history, given in Figures |6.3| and |6.4[ respectively, indicates 
similar behaviour to the three-dimensional case. The choice in the shape parameter 



is as in (5.1 1 with K = 3. We remark on the favourable complexity of the MLSKI al- 
gorithm compared to RBF/MLRBF as d grows. Indeed, the MLSKI algorithm is able 
to compute highly accurate interpolants, as large N can be achieved using moderate 
computational time. 

6.4. Experiment 4. It is easy to see that d-dimensional Gaussian kernels are 
tensor-products of one-dimensional ones. Hence, it is interesting to also consider SKI 
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and MLSKI with non-Gaussian kernels. 

To this end, we consider the problem of interpolation of uf3.d using Wendla nd's 
compactly the supported kernel (j>3,2(r) ■— (1 — cr) f j_(35c 2 r 2 + 18cr + 3). In Figure 



6.5 



the convergence history of the various interpolation methods based on the compactly 



supported kernel 4>3,2 are given. The choice in the shape parameter is as in (5.1| 
with K = 3, resulting in safe conditions numbers for all methods. Interestingly, the 
SKI method seems to converge very slowly, while the MLSKI algorithm appears to 
perform well. 

Finally, we consider the problem of interpolation of ufad and u qua d using the 
inverse multiquadric kernel (Pimq{t) := (1 + c 2 r 2 ) -1 / 2 with the above choice of the 
shape parameter. The results are given in Figures |6.6| and |6.7[ respectively. Again, 
for the case of u qua d the SKI methods seems to be performing poorly, compared to 
the fast convergence of MLSKI. 

7. Concluding remarks. A multilevel kernel-based interpolation method, suit- 
able for moderately high-dimensional interpolation problems on carefully structured 
grids has been proposed. The key idea is the use of hierarchical decomposition of the 
data sites with anisotropic radial basis functions used on each level, solving a num- 
ber of smaller independent interpolation problems. MLSKI appears to be generally 
superior over classical radial basis function methods in terms of complexity, run time 
and convergence, at least for large data sets when d = 3, 4. It is expected that good 
convergence can be obtained for d = 5 also. 

Currently, the choice of data sets (sparse grids) is highly structured. Some pre- 
liminary numerical experiments for SKI on mildly perturbed sparse grids, presented 
in [JB], indicate that the SKI method converges with the same rate, albeit with a 
somewhat larger constant. Perhaps a more robust methodology for extending the 
applicability of SKI/MLSKI methods to scattered data is the pre-computation of the 
values on the corresponding sparse grid data sites via local interpolation. The exten- 
sion of the SKI/MLSKI method to more general geometries could possibly be handled 
either by introducing fictitious gridded data sites with suitable data values, taking into 
account the nature of the data, or by conformally mapping the computational domain 
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Fig. 6.6: Convergence of 4>imq for uf^d- Error evaluated at 194,481 Halton points. 
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Fig. 6.7: Convergence of 4>imq for u q 
Halton points. 




■ RBF 
MLRBF 

■ SKI 

- MLSKI 



for d 



Time 

Error evaluated at 194,481 



The discussion in this work has been confined to strictly positive definite kernels. 
The implementation of MLSKI with conditionally positive definite kernels is subject 
to ongoing work; some preliminary results can be found in [46] . 
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